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Abstract 

Numerical algorithms are proposed for simulating the Brownian dynamics of charged particles in an 
external magnetic field, taking into account the Brownian motion of charged particles, damping effect and 
the effect of magnetic field self-consistently. Performance of these algorithms is tested in terms of their 
accuracy and long-time stability by using a three-dimensional Brownian oscillator model with constant 
magnetic field. Step-by-step recipes for implementing these algorithms are given in detail. It is expected 
that these algorithms can be directly used to study particle dynamics in various dispersed systems in the 
presence of a magnetic field, including polymer solutions, colloidal suspensions and, particularly complex 
(dusty) plasmas. The proposed algorithms can also be used as thermostat in the usual molecular- dynamics 
simulation in the presence of magnetic field. 
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I. INTRODUCTION 



Brownian Dynamics (BD) simulation method for many -body systems of particles immersed in 



a liquid, gaseous or plasma medium 



can be regarded as a generalization of 



the usual Molecular Dynamics (MD) method for many-body systems in free space. While the MD 
method is based on Newton's equations of motion^he BD method is based on their generalization 
in the form of Langevin equation and its integral ||8|] : 
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-F + A{t), 
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where, as usual, m, v and r are, respectively, the mass, velocity and position of a Brownian par- 
ticle, whereas F is a systematic (deterministic) force coming from external fields and/or from 
inter-particle interactions within the system. What is different from Newton's equations is the ap- 
pearances of dynamical friction, —7V, and random, or Brownian acceleration, A(t). These two 
force components represent two complementing effects of a single, sub-scale phenomenon: nu- 
merous, frequent collisions of the Brownian particle with molecules in the surrounding medium. 
While the friction represents an average effect of these collisions, the random acceleration repre- 
sents fluctuations due to discreteness of collisions with molecules, and is generally assumed to be 
a delta-correlated Gaussian white noise. The friction and random acceleration are related through 
a fluctuation-dissipation theorem which includes the ambient temperature, therefore guaranteeing 
that a Brownian particle can ultimately reach thermal equilibrium within the medium ||8|]. 

The Langevin equations, Eq. Q), can be numerically integrated in a manner similar to New- 
ton's equations in the MD simulation, which gives rise to several algorithms for performing BD 
simulation, such as the Euler-like [1], Beeman-like [2, 3]> Verlet-like [3], and Gear-like Predictor- 



Corrector (PC) methods [11511 . as well as a wide class of Runge-Kutta-like algorithms (see, e.g.. 
All those methods were used successfully to stu dy p roblems in various dispersed sys- 



tems, sue 



1 as po 



plasmas jll . 
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ymer solutions [|9D, colloidal suspensions HlOll and, in particular, complex (dusty) 



Recently, there has been a growing interest in studying the dynamics of dust particles and dust 



clouds in both unmagnetized jln . 



1811 and magnetized dusty 119 . 



20, 



21 



,12211 plasmas. The topics 



studied so far include, besides Brownian motion of a dust particle in an unmagnetized plasma 



nlm . also the gyromotion of a single dust particle and rotations [19] of dust clouds M 



21, 
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and clusters [1231 



in a magnetized plasma. There have also been several theoretical 
proposals for studying waves and collective d yna mics [ 27I1 in such magnetized plasma systems, in 



which dust particles are fully magnetized [|28 



291 



30, 



31 



32ll . However, exploring those proposals 



in the laboratory does not seem to be quite feasible as yet, due to many constrains [|28 



29 



30 



320. Therefore, it is desirable to have algorithms for numerical experiments that can validate the 
existing theories, on one hand, and that can serve as a guide for future laboratory experiments, on 
the other hand. Since there are no such algorithms, to the best of our knowledge, we propose here 
a few new BD algorithms, which treat an external magnetic field in the simulation in a manner 
consistent with the Langevin dynamics, Eq. ©. 

The manuscript is organized in the following fashion. In Sec. UIl we present the general for- 
mula for a BD simulation with a constant magnetic field. Detailed implementations to the Euler-, 
Beeman- and Gear-like methods are given, respectively, in Sections Unl |IVl and |Vl Concluding 
remarks are contained in Sec. I VIII 



II. GENERAL FORMULA FOR BD SIMULATIONS 



The dynamics of a charged Brownian particle in a constant external magnetic field B is de- 
scribed by introducing the Lorentz force in the Langevin equation [|7i ISD 
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7V + — F + -^vxB + A(t), 



mc 



(2) 



where Q is the charge on the particle and c is the speed of light in vacuum. 

As usual, certain assumptions must be made about the deterministic force F in order to con- 
struct a meaningful algorithm for a many-particle simulation from the above equation. A common 
approach [1, 2, 2, is to assume that F is only an explicit function of time t. Thus, a Taylor series 
of F or, equivalently the deterministic acceleration, a = F/m, can be written as 



a(t) = a(0) + a(0)t + ^8L{0)t^ + ^a(0)t^ + 



(0)f 



(3) 



where a*^") represents the nth-order time derivative of a. There are various ways to derive formulas 
for conducting a BD simulation from the Langevin equation, Eq. ([2]), based on the above Taylor 



series. We shall adopt the strategy outlined in Refs. 

mm 

, and more recently implemented 
in Refs. [|7i [SD, as it is simple and straightforward, especially for readers with some simulation 
background but without much background, or interest in stochastic calculus. 

The Langevin equations, Eq. Q, may be integrated analytically in a short time, based on an 
adopted truncation rule for the series in Eq. ([3]), thereby giving an updating formula for a BD 
simulation. (We note that detailed technique for integrating the Langevin equation is particularly 
well described in Refs. fllsl].) The resultant formulas are actually expressions for the two random 
variables, v(t) and r(t), which, under the assumptions that the Brownian acceleration in Eq. Q is 
a Gaussian white noise and that 7 is constant, turn out to be normally distributed random variables, 
according to the normal linear transform theorem [ItI,^]. Therefore, the Cartesian components of 
the velocity v = {v^, Vy, Vz} and position r = {x, y, z} vectors for a Brownian particle can be 
expressed in terms of their respective means and variances 



v^{t) = mean{t;,(t)} + Vvar{t;,(t)}iV^(0,l), 

a{t) = mean{a(t)} + v/var{a(t)} iV^(0, 1), (4) 

where a takes values x, y and z, and A^(0, 1) is a shorthand notation for the standard normal ran- 
dom variable having zero mean and unit variance, the so-called unit normal 13,181]. The superscripts 
attached to the components of random vectors N"' = {N^,N^, NJ} and N'" = {N^, N^, N^}, 
appearing in Eq. dU, indicate that those two sets of unit normals are associated, respectively, with 
the velocity and position of the Brownian particle. We note that the Cartesian components of the 
vector N'^ are mutually independent, as are the components of the vector N"", but the vectors N"^ 
and N"" are correlated, as will be shown below. 

We should emphasize the importance of Eq. (H)) because it provides a general updating formula 
for the BD simulation and also enables substantial simplifications in the subsequent design of our 
algorithm. Now, solving the stochastic differential equation, Eq. Q, and obtaining the two random 
variables v(t) and r(t), is simply reduced to determining the two sets of deterministic quantities, 
the means and variances of v(t) and r(t), as well as their covariances, assuming an appropriate 
truncation scheme for the deterministic acceleration, Eq. ([3]). To simplify the notation, in the 
following we shall denote the velocity and position means, respectively, by (va) = mean{va} and 
(a) = mean{a}, and we shall use standard deviations a instead of variances. 
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A. Variances and covariances 



We begin with variances and covariances because they do not depend on the form of determin- 
istic acceleration, a(t). 

It should be noted that random increments in the velocity do not depend explicitly on the 
magnetic field and are therefore isotropic. The corresponding variances are then given by 

\ar{iKj.} = \ar{vy} = varjf^} = var{t>} = cr^, where 



^. = \/^(l-e-2^0, (5) 
V m 

with ks being the Boltzmann constant and T the temperature of the medium. [It should be noted 
that in an equilibrium between the Brownian particle and the medium, T is also the kinetic tem- 
perature of the Brownian particle. However, Brownian particles may have kinetic temperature that 
is different from T, which opens a possibility of simulating non-equilibrium processes by using 
the BD]. 

However, random displacements of the position do depend explicitly on the magnetic field 
and therefore are non-isotropic. Let us assume B = {0, 0, B} and define Q = QB/{cm) to be 
the gyrofrequency of a Brownian charged dust particle. Then, we have var{a;} = var{y} = aj_ 
and var{^} = a^, where a± and a\\ are, respectively, standard deviations of the position in the 
directions perpendicular and parallel to the external magnetic field 13, l8|] . In particular, we find 
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Furthermore, we have (6 x 5)/2 covariances 
covariances are co\{vx, x}, co\{vy, y}, co\{vx, y}, co\{vy, x} and coy{vz, z} 
the definition cov{X, Y} = (XY) — {X) (Y) for the covariance of two random variables, X and 
Y). In particular, we obtain lAm 



most of which are zeros. The non-zero 
flS] (here, we use 
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K = cow{x,v^} = cov{y,Vy} = t—\ S' (l + e"^^* _ 2e"^* cos fit) , 

m '-ft'-f'^ + \ r ' 

H = covlv, vA = — cov|x, Vy} = t — — ( 1 — e"^'^* — 2e~^^— sin fit ) , 

^ I ' 2/i m 7^72 + fi2 V fi / ' 

kT 

L =cov{v,,z} = t^(l-2e-"" + e-^""). (7) 
m'jt ^ ' 

We note that, when i? ^ 0, i.e., fi ^ 0, one has H ^ and K L, so that a± a y, recovering 

the results for systems without magnetic field [jsl]. 



B. Mean values 



Given the expression for a(t) in terms its Taylor series in Eq. dB]), one can obtain the means 
(v) = {(^'x), (vy), {vz)} and (r) = {(x), {y), {z)} in a number of ways by integrating Eq. Q. For 
simplicity, we follow here Lemons' methodology DTI, ISD, in which Eq. ^ is reduced to a set of 
deterministic (ordinary) differential equations by taking expectation values of both sides 



d(v) 
dt 

dt 



-7(v)+a(t) + -^(v)xB, 
mc 



(8) 



(v), 



and using the fact that the Brownian acceleration A{t) is a Gaussian white noise with mean 0. 

Equation ([8]) can be integrated analytically by assuming that initial conditions fq, vq, ao, ao, 
ao, a'o, . . . , SL^\ are known at t = 0. One ends up with 



ro + Ii ■ vot + I2 ■ + I3 ■ ko^ 

j-2 



(r) 

(v) = lo ■ vo + Ii ■ aot + I2 • aot^ + I3 ■ kot 
With the assumption B = {0, 0, 5}, the matrices I„ are given by 
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for n = 0, and 
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for n > 1. 

It is interesting to note that the coefficients satisfy very simple recursive relations 



(12) 



(M") = &„-it , 

dt 



(dnf) = dn-lt 



n-1 



(13) 



A special yet important case would be that of = 0, i.e. = 0, and consequently 6„ = and 
Cn = dn, so that In bccomes a diagonal matrix with the diagonal elements being dn and, Eq. ^ 
reads as [|l5l] 



(r)B=o = ro + rfivot + rfaaot^ + rfgaot^ + ■ ■ ■ + c^naS," ^^t" + ■ ■ 
(v)b=o = (iovo + rfiaot + (i2aot^ + dsaot^ H h ci„,a[,""^^t'' + 



(14) 



Further, when 7 — * 0, c„ and dn reduce to the usual coefficients of a Taylor series. 



C. General updating formula 



With the variances, covariances and means known, the updating formula for a BD simulation 
can be further written in a compact form as 13, [SO 



v(t) = (v) + (T,Ni(0,l), 

r(t) = (r)+I-Ni(0,l) + J-N2(0,l^ 



(15) 
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V^" 

The unit normal vectors Ni(0, 1) and N2(0, 1), besides each having statistically independent 
Cartesian components, are also by design statistically independent from each other, and may be 
generated by the Box-MuUer method [330 on a computer. 

As expected, when B = 0, I and J become diagonal matrices with their diagonal elements 



being L/a^, and y'aj| — L'^/cr^, respectively, in which case Eqs. (fT5l) become 

VB=o(t) = (v) + a,Ni(0,l) 



rB=oW = (r) + — Ni(0,l 




^N2(0,l). 



(18) 



It is worthwhile to point out that the updating formulae, Eqs. (1151) and (1181) . can be actually 
interpreted as a two-step algorithm: first, one calculates the means of the velocity and position at 
time t, and secondly, one adds explicitly random displacements (ERDs) of velocity and position. 
This is very important in the sense that the first step is nothing more than solving deterministic 
Newton's equations with damping, so in principle many algorithms suitable for the MD simulation 



3411 ) can be used here. On 



(for example, the Verlet, Beeman and even multi-step PC algorithms [|6L 
the other hand, the second step is independent of the first one, and can always be performed at the 
end of the time step as a correction to the previous step. 

However, Eqs. ^ and (fT4l) are still not very suitable for computer simulations, because they 
involve high-order derivatives of acceleration, which are not directly accessible in simulations. 
Further assumptions must be made about a(t), and different ways of handling that issue then 
result in different simulation methods, such as the Euler-like, Beeman-like, Verlet-like and Gear- 
like methods. 



III. EULER-LIKE METHOD 



In the Euler-like method yj], it is assumed that F = F(0), or a = a(0) = constant in the time 
interval [0, t]. Under this assumption, the velocity and position means are, respectively, given by 

{r)EL = ro + Ii ■ vot + h ■ aot^ 

(v)ijL = lo ■ vo + Ii ■ aot, (19) 



8 



where the subscript EL in Eq. (fT9l ) is added to indicate reference to the Euler-like method. The 
special case of zero magnetic field can be simplified to yj] 

(r)sL,B=o = ro + (iivot + c?2aot^ 

{^)el,b=o = dovo + diaot. (20) 

The procedure for implementing the simulation by using Euler-like method is as follows: 

(a) Evaluate Eq. (fT9l ) [or Eq. (|20l ) in the case of zero magnetic field] for the means of velocity and 
displacement, given the initial conditions vq, fq and ag at t = 0. 

(b) Generate two independent random vectors, Ni(0, 1) and N2(0, 1), which follow the standard 



normal distribution (this may be usually done by using the Box-MuUer method 113311 ). and sub- 
stitute them into Eq. (fTSl) [or Eq. ([TSl) in the case of zero magnetic field], along with the means 
obtained in the first step. 

(c) After the above two steps, a force evaluation is done, which will provide an initial condition 
for ao, along with the updated position and velocity, to be used in the next step. 



IV. THE BEEMAN-LIKE METHOD 



In the Beeman-like method ||2, 



i 



one has 



{r)BL = ro + Ii ■ vot + h ■ aot^ + I3 • ao^^ 

(v)bl = lo ■ vo + Ii ■ aot + I2 • aot^ + I3 ■ aot^ (21) 

Note that, according to the Beeman algorithm |2, ^, terms through first order in t, i.e., a(t) = 
ao + agt, are kept in the expression for (r) bl, and terms through second order in t, i.e., a(t) = 
ao + act + Ifaot^, are kept in (v)bl- 



However, the above expressions are not yet suitable for a one-step simulation method [1340, as 



they contain derivatives of the acceleration, which have to be replaced by a finite difference for- 
mula in terms of a(0), a(— t) and a(t). Thus, we obtain the schemes of the Beeman-like algorithm 



{r)BL = ro + la ■ vot + If, ■ aot^ + Ic ■ a_ft^, 
{v)bl = Id • Vo + le • att + 1/ • aot + I^ ■ a_it, 



(22) 
(23) 
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where at = a(t), a t = a(— t), and 



Ia = Il, Ife = l2+l3, Ic = -l3 

h = lo, le = I2 - I0I3/I1, 1/ = Il - I2 + 2I0I3/I1, Ig = -hh/h- 

[Note that, here, matrix operations should be understood as direct element operations, e.g., 

(Iol3)ji = lo.ijls.ii-] 

In the special case of zero magnetic field, Eqs. (123]) and (|22l) can be simplified to ||2|] 

(r)BL,B=o = r(0) + c,v(0)t + ci,a{0)f + c,a{-t)t^, (24) 
(v)i3L,B=o = Crfv(O) + Cea(t)t + c/a(0)t + Cga(-t)t, (25) 

with 

Ca = C?i, Cb = d2 + rfs, Cc = -rfs 

Cd = do, Ce = d2- dod-i/di, c/ = cii - 4 + 2c?oC^3/c^i, Cg = -d^d^/di. 

Implementation of the Beeman-like method is also based on Eq. ( fTSl ). and the detailed simula- 
tion procedure is listed in the following: 

(a) Evaluate Eq. (|22|) [or Eq.(l24l)1 for the mean of displacement, given the initial conditions vq, tq 
and ao at t = 0, and a(— t) at time —t. 

(b) Evaluate a new deterministic acceleration a{t) based on the new position obtained in the pre- 
vious step. 

(c) Evaluate Eq. (|23] ) [or Eq.(|25l)1 for the mean of velocity, given the initial conditions vq and ao at 
t = 0, and a{—t) at time —t, as well as the newly obtained acceleration a(t) in the previous step. 
[Note that, if one uses a periodic boundary condition, it should be applied before this step.] 

(d) Generate two independent random vectors, Ni(0, 1) and N2(0, 1), which follow the standard 
normal distribution, and substitute them into Eq. ([T5l) along with the means obtained in the steps 
(a) and (c). 

(e) Apply periodic boundary condition again if using it. 

As one can see, the difference compared to the Euler-like method, is that here one has to do 
the force evaluation after the position update, but before the velocity update, in every time step, 
and one has to store the force value of the last step to be used for a(— t). In addition, one needs to 
apply the periodic boundary condition twice. 
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Let us note that the Verlet-like method yD was based on similar assumptions for a(t) as the 
Beeman-like method, and it was proven y20 that the Verlet-like method is numerically equivalent 
to the Beeman-like method in the position, while the latter seems to have a better accuracy in the 
velocity. Therefore, extensions of the Verlet-like method are not presented here. 



V. GEAR-LIKE PREDICTOR-CORRECTOR METHOD 



A Gear- like Predictor-Corrector (PC) method for BD [llSll can be constructed in a direct analogy 



with the Gear method for MD simulations. Our Gear-like method also includes three stages. 



340, as in the MD simulation, but the 



namely, predicting, force evaluating, and correcting 
difference here is that one has to add explicit random displacements of velocity and position by 
using Eq. (fT5l) at the end of time step to complete the BD simulation. The basic procedure goes as 
follows. 



A. Predicting 

In the predicting stage, one has 

(r)^ = ro + Ii ■ vot + h ■ aot^ + I3 ■ aot^ + I4 ■ aot^ + I5 ■ aot^ 
(v)^ = lo ■ vo + Ii ■ aot + I2 ■ aot^ + I3 ■ aot^ + I4 ■ ao^^ 
= ao + aot + ^'^ot^ + ^'^'0^^' 

a^ = 3.0 + aot + — a'ot^, 

a^ = 3.0 + a'ot, 

a^ = ao, (26) 

where the superscript P indicates that these are quantities in the predicting stage. For simplicity, 
we have dropped in the above all derivatives of a(t) higher than the third order, but we note that 
extensions to higher orders are quite straightforward. One notices in Eq. (|26|) that we have used 
Eq. ^ for the means of position and velocity, instead of using Taylor series for position and 
velocity which usually provide a basis for implementing the Gear method in the MD simulation 



34|] . In the case of zero magnetic field, one simply needs to replace the means of velocity and 
position given by Eq. (fT4l) . instead of Eq. The rest of the algorithm (derivatives of the force) 
is essentially the same as in the MD. 
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B. Force evaluating 



In the next step, the predicted position r^(t) is used to obtain a new force, or deterministic ac- 
celeration a(t), and a difference between the predicted acceleration a^{t) and the new acceleration 
a(t) is calculated as 



Aa = a(t) - a^(t). 



(27) 



It can be seen that this step is exactly the same as the one normally used for the Gear method in 
MD l|6,34]. 



C. Correcting 

In the correcting stage, the above difference term is further used to correct all predicted posi- 
tions and their "derivatives", thus giving 

{rf = {rf + 2aol2 ■ AR, 
(v)^t = (v)^t + aili ■ AR, 
+ aaAR, 

+ ttgAR, 

+ a4AR, 

+ asAR, (28) 
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where 

AR = (29) 
and the coefficients are given in the Table HI Note that the table is simply a reproduction of 



those appearing in Refs. [[a [SJ], and is given here for completeness. By using parameters in 



different columns of that table, one can achieve 3rd-, 4th-, and 5th-order (or 4-, 5- and 6-value) 
Ifi, 34] Gear-like algorithms for the BD simulation. Note that the first two equations in Eq. (|28]) 
are slightly different from those in the MD in order to restore the damping effect on deterministic 
acceleration, and to maintain consistence with the corresponding terms in Eq. Q [or the first two 
equations in Eq. (|26|) 1 as well. In the case of zero magnetic field, one simply replaces Iq and Ii in 
Eq. (l28l) by do and di, respectively. 
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TABLE I: Coefficients used in correcting stage of Gear-like PC method. Note that the table is simply a 
reproduction of those appearing in Refs. jq. 13411 . 



ai 3th-order 4th-order 5th-order 

ao 1/6 19/120 3/16 

ai 5/6 3/4 251/360 

02 1 1 1 

as 1/3 1/2 11/18 

04 1/12 1/6 

05 1/60 



D. Adding explicitly random displacements (ERDs) 

To complete the BD simulation, we have to use the updating formula, Eq. (fT5l) . to add ERDs of 
the velocity and position. It should be noted that now the corrected value (r)*^ and (v)*" must be 
used, respectively, in the places of (r) and (v) in Eq. (fT5l) . 

We summarize the basic steps for implementing the Gear-like PC method: 
(a) Eq. (|26l ) is used to calculate predicted values of the position, velocity, acceleration and its 
derivatives, given the initial conditions, tq, vq, ao, ao, ao and a'o, at t = 0. Note that, for the very 
first few steps of the simulation, ao, ao and a'o are undefined. The simplest way to get around 
this issue is to simply set all of them to zero at the very first step, and their values then will be 
updated during subsequent iterations. A better way would be to start the simulation by using a 



Runge-Kutta procedure for the first few steps ||34|]. However, neither of these alternatives will 



have any significant effects on the results in real many-particle simulations. 

(b) Evaluate new acceleration a{t) by using the predicted position (r)^, and calculate its differ- 
ence with the predicted value a^(t) by using Eq. (ITTI) . [Note that, if one uses periodic boundary 
condition, it should be applied before the force evaluation.] 

(c) Correct the predicted values of the position, velocity, acceleration and its derivatives by using 
Eq. dH. 

(d) Generate two independent random vectors, Ni(0, 1) and N2(0, 1), which follow the standard 
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nomial distribution, and substitute them into Eq. (fTSi) together with the corrected values (r)*" and 

(e) Apply periodic boundary condition again if using it. 

The above is a basic procedure for using the Gear-like PC method for the BD. One might have 
noticed that the formulas, as well as the simulation procedure, are quite similar to those used in 



the Gear method for the MD based on Newton's equations [I6|,l34|], apart from our use of Eqs. (I26|) 
and (|29l) to express the velocity and position, as well as the addition of ERDs at the end of every 
time step. Also, when B = and 7 — 0, the Gear-like method goes over to the Gear method for 
MD simulation. 

VI. TESTING THE ALGORITHMS 



In the above we have developed numerical algorithms for simulating Brownian dynamics of 
charged particles in an external magnetic field. But how accurate are they in describing the actual 
Brownian motion? To answer this question, we present in this section some simple computational 
examples as testing cases and compare the performances of different algorithms presented above. 
For simplicity, we shall occasionally denote the Euler-like, Beeman-like, Verlet-like and Gear-like 
methods by EL, BL, VL, and GL, respectively. 

As the simplest test cases, one could adopt comparison of the numerical results with the ana- 
lytical results for some simple model problems, such as the classical harmonic oscillator, which 



had been used extensively to test algorithms for the MD simulation (see for example Ref. ll34ll 
and ^Bsj] for nice reviews). We follow here the same logic and employ the model of a three- 



dimensional (3D) Brownian-harmonic-oscillator (BHO) in an external magnetic field ll40ll . for 
which the Langevin equation (O is reduced to 

dv 

—M. = -^yy-ujly-nv^ + Ay{t), 
du 

-£ = -^y^-ujlz + A,{t). (30) 

Since the magnetic field is in the z direction of the Cartesian coordinate system, two independent 
processes take place. In the z direction, the magnetic field does not have any effect on Brownian 
motion, so the Brownian particle behaves like a one-dimensional stochastically-damped harmonic- 



oscillator (Si, 



36h . while in the directions perpendicular to the magnetic field, i.e., in the xy plane, 
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FIG. 1 : (Color online) Mean position (left) and mean velocity (right) of a Brownian oscillator in magnetic 
field during 200 time units, with 7 = 0.02, ujq = l/\/2, ft = 0.5 and At = 0.02. Initial conditions are 



xq = 1.0, yo = 0, zq = 1.0 and v^o = Vyo = Vzo = 0. Solid lines are the result of analytical solutions 14011 . 
and circles are the numerical results calculated by using the GL-5 method. 

Brownian motion is much more complicated because of the coupling between the motions in the 
X and y directions via magnetic field. It is curious to note that, although the analytical solution 
for a BHO without magnetic field had been known for well over a half of the century jBg] , the 
problem of BHO in an external magnetic field has been solved analytically only very recently by 
Jimenez- Aquino et a/. ll40ll . 

Of course, the dynamics a BHO in magnetic field can be also traced by using the above de- 
scribed numerical methods. It is important to realize that, although the physical model for a BHO 
in magnetic field, Eq. (|30l ), involves a very simple deterministic acceleration which depends on 
the particle's position, it nevertheless provides a good test for our assumption that this acceleration 
depends on time only, and to examine the effects of various truncation schemes for the Taylor 
series representation of this acceleration, i. e., Eq. ([3]). The results of our numerical simulations of 
the BHO in a magnetic filed will be compared with the corresponding explicit analytical solutions, 
for which we refer the reader to the original reference ||40|]. By doing so, one can evaluate the 
accuracy and performance of the numerical methods and validate the assumptions made in their 
derivation. The results will serve as a basic reference for future simulations of more complicated 
systems, such as magnetized dusty plasmas 113 511 . 
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FIG. 2: (Color online) Components of the mean position (left) and mean velocity (right) in the xy plane of 
a Brownian oscillator in magnetic field during 200 time units, with 7 = 0.02, ojq = l/\/2, fl = 0.5 and 
At = 0.02. Initial conditions are the same as in fig.[I] Solid lines are the result of analytical solutions 14011 . 
and circles are the numerical results calculated by using the GL-5 method. Dashed lines are projections of 
the analytical results on the xy plane. 



According to the above discussion, and particularly referring to Eq. (fTSl ), the task of a BD 
simulation is simply to predict the position r(t) and velocity v(t) of a Brownian particle at time 
t, given the set of initial conditions at time 0. Since r(t) and v(t) can be obtained numerically 
in two-steps, first, by calculating the means of the velocity and position at time t and, second, by 
adding the ERDs of velocity and position, the performance of a BD simulation will be examined 
by testing the accuracy of both steps. We begin by testing the first step, i.e., calculating the means. 



A. Inaccuracy in calculating the means 

Without any loss of generality, here and in the subsequent simulations we set fc^T = 1, m = 1 
and ujq = and choose the initial conditions to be = 1-0, yo = 0, zq = 1.0 and v^o = 

VyO = VzO = 0. 

We first present in Fig. [T] examples of full 3D trajectories of the mean position and velocity for 
a BHO in 200 time units, with 7 = 0.02, Q = 0.5 and the time step size At = 0.02. Solid lines are 



the results of the analytical solutions [Eqs. (B23)-(B34) in ll40n i. while the circles are numerical 
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FIG. 3: (Color online) Deviations of the mean x (solid lines) and mean Vx (dashed lines) values from the 



analytical solutions 14011 using different methods (from top to bottom: GL-5, GL-4, GL-3, BL and EL), with 
7 = 0.02, u>o = l/\/2, ^ = 0.5 and At = 0.02. Initial conditions are the same as in Fig. [T] Note the 
magnitudes of the errors. 



results calculated by using the GL-5 method. One might be particularly interested in the motion in 
the directions perpendicular to the magnetic field, which is shown in Fig. [2] as the time evolutions 
of the mean position and the mean velocity in the xy plane. We have found that details in the 
trajectory patterns strongly depend on the initial conditions and on the values of ujq and il, but the 
general tendency in trajectories is the same. Without Brownian acceleration, the initial energy of 
the oscillator would be sooner or later consumed by the damping, and the oscillator would come 
to rest at X = y = z = after long enough time. One sees from Figs. [T] and [21 that the numerical 
results agree very well with the analytical solutions. Good visual agreements were also found with 
numerical results of the GL-4, GL-3 and BL. However, quantitatively, they are quite different, as 
is shown in the following. 
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FIG. 4: (Color online) Deviations of the mean z (solid lines) and mean Vz (dashed lines) from the analytical 



solutions 14011 using different methods (from top to bottom: GL-5, GL-4, GL-3, BL and EL) in 200 time 
units, with 7 = 0.02, ljo = l/\/2, = 0.5 and At = 0.02. Initial conditions are the same as in Fig. [T] 
Note the magnitude of the errors. 

Figures [3] and |4] display deviations (differences) of the numerical results for the means of the 
position and velocity from the corresponding analytical results ll40ll in 200 time units, under the 
same conditions as in Fig.d] Only deviations in the x and z directions are shown in Figs. [3] and IH 
respectively, as those in the y direction are essentially the same as those for x, apart from a phase 
shift. The oscillatory patterns of deviations in the position and velocity approximately resemble 
those of the full solutions in Fig. [H but they have much smaller amplitudes. One can observe 
the differences in magnitude of deviations for different methods, with the GL-5 method having 
the smallest deviations and therefore highest accuracy in both the x and z directions, while the 
EL method exhibits the largest deviations, as can be expected from our previous comparisons of 
these methods. It should also be noted that the performance of these methods is different in the x 
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FIG. 5: (Color online) The largest deviation in the position (x-component) from the corresponding ana- 
lytical solution 1 4^ in the first 20 time units versus il, showing the dependence of inaccuracy of different 
methods on the magnetic field, with = l/\/2, 7 = 0.02 and At = 0.02. 



and z directions. One sees generally higher accuracy in the z direction for the Gear-like methods, 
while the BL and EL methods have similar accuracies in the two directions. This indicates that the 
presence of the magnetic field also affects the accuracy of the computation. 

A more detailed analysis of the dependence of the accuracy on the magnitude of magnetic 
field is shown in Fig. [51 where the largest deviations of the position in the x direction from the 
corresponding analytical solution! 40|] are recorded in the first 20 time units (which is a convention 



for the measure of accuracy DalMl]) and are plotted versus for different methods. One observes 
that, with the increase of the magnetic field intensity, the largest deviation of the EL method 
remains almost constant, while that of the BL method initially stays constant, but drops slightly 
when r2 > 1.0. The tendency observed in the Gear-like methods is the opposite. The initially 
excellent accuracy deteriorates with increasing magnetic field. The one with the highest accuracy, 
i.e., the GL-5 method is affected the most, as is shown. A similar scaling rule applies to simulations 
with other time steps and damping rates, as discussed next. 

We make comparisons involving different time steps because that is always a key issue in both 
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FIG. 6: (Col 



or online) The largest deviation in the position (z-component) from the corresponding analytical 



solution 14011 in the first 20 time units versus time step size At, with = ^/V^ and = 0, for different 



methods and different 7 values. 



the MD and BD simulations. Fig. [6] shows the largest deviations (defined by the largest deviation 
in the first 20 units) of the position in the z direction versus the size At of the time step for different 
methods and different friction coefficients 7. This figure is plotted in a double logarithmic scale, 
and the curves are nearly straight lines. The slope of these lines is called the apparent order ll34ll . 
illustrating the dependence of the error on the time step size. Namely, if the error is found to be 
proportional to At^, then the exponent p is the apparent order. It is found that for very small 7, for 
example 7 = 0.01, as shown in Fig.[6l the apparent orders are pel ~ 1, Pbl ~ 2, pgl-s ~ 3.5, 
~ 4.2 and Pgl~5 ~ 4.6, respectively for the Euler-like, Beeman-like, 3rd-, 4th-, and 5th- 
order Gear-like methods. These values are very close to those from the MD simulations where 
damping is absent ll34ll . which proves the consistency of our computation. When 7 increases, 
the absolute value of the error for Euler-like and Beeman-like methods decreases, while pel and 
Pbl remain almost unchanged. On the other hand, Pgl-s, Pgl-a and pcL-b slightly decrease 
with increasing 7. For example, for 7 = 1, as shown in Fig. [6l the Gear- like methods have the 
worst performance in accuracy: their apparent orders become now Pgl~3 ~ 3.1, Pgl-a ~ 3.9 
and Pgl-5 ~ 4.0), but the magnitudes of their errors are still much smaller than those of the 
Beeman-like and Euler-like methods. 

All the above tests show that the numerical methods, especially the Gear-like methods and 
Beeman-like method, can describe the mean values of the movement of a BHO with sufficiently 
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FIG. 7: (Color online) Position (left) and velocity (right) of a Brownian oscillator in a magnetic field during 
200 time units, with 7 = 0.02, ljq = l/\/2, 17 = 0.5 and At = 0.02. Initial conditions are the same as in 
Fig.[T] Solid lines are the result of analytical solutions |4ot, and circles are the numerical results calculated 
by using the GL-5 method. 



high accuracy. 



B. Total inaccuracy 

The tests carried out in this sub-section are similar to those in the previous one, but with the 
addition of the ERDs in both the position and velocity. We first show in Fig. |7]full trajectories of 
the BHO in 200 time units, under the same condition as in Fig.[TJ Again, solid lines show results of 



the analytical solutions 114011 . while circles are the numerical results calculated by using the GL-5 
method. Motion of the BHO in the xy plane under the influence of magnetic field is shown in Fig. 
[8l One observes in both figures that the numerical results of the GL-5 method again agree very 
well with the analytical results. 

Before giving a more quantitative analysis of the full deviations including the ERDs by using 
different methods, we comment on the accuracy in variances and covariances, Eqs. © and ([7]), 
which will help to better understand the deviations in full trajectories. In deriving Eqs. © and 
(|7]), we used the assumption that the deterministic force is an explicit function of time only, and 
we truncated the Taylor series representation for that force. However, in many realistic problems 
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FIG. 8: (Color online) The xy components of the position (left) and velocity (right) of a Brownian oscillator 
in a magnetic field during 200 time units, with 7 = 0.02, ujq = l/\/2, = 0.5 and Ai = 0.02. Initial 
conditions are the same as in Fig. [T] Solid lines are the result of analytical solutions ll40ll . and circles are 
the numerical results calculated by using the GL-5 method. Dashed lines are projections of the analytical 
results on the xy plane. 

and in the models such as the BHO, the deterministic force depends explicitly on the particle 
position only. While we have seen in the previous sub-section that this assumption can provide a 
sufficiently accurate description of the means of the position and velocity for a BHO, a question 
remains as to how does this assumption affect the calculation of variances and covariances, i.e. the 
explicitly random part of displacements. 

This is addressed in Fig. |9] which depicts relative deviations of the variances and covariances 
computed by the GL-5 method from the corresponding analytical results [Eqs. (26)-(34) in [|40ll 1 
versus the time step size At for 7 = 0.02 and Vt = 0.5. (Note that, in the specific case of zero 



magnetic field, the analytical results of [14011 are identical to those of Refs. |18|,|36|].) All deviations 
are seen to increase in Fig. |9] with increasing At. In the log-log scale, all results form a cluster of 
parallel straight lines with a slope of about 2.5, indicating an apparent order of approximately 2.5. 
Judging by both the magnitude of deviations and by the apparent order, the accuracy for variances 
and covariances displayed in Fig.|9]is much lower than the accuracy for the means shown in Fig. 
|6l One would expect that the total accuracy of an algorithm will be largely determined by the 
accuracy of that part of the algorithm which has the lowest accuracy. 
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FIG. 9: (Color online) Deviations of variances and covaiiances from the conesponding analytical solutions 



1 4011 versus the time step size At, for u>q = l/\/2, 7 = 0.02 and Q. = 0.5. 



As before, we next carry out simulations over certain periods and record full deviations (means 



plus ERDs) of the position and velocity from the corresponding analytical results 114011 . All results 
are assembled in Fig. [10] for different methods in 200 time units, with At = 0.02, 7 = 0.02 
and f2 = 0.5. One sees that, for all Gear-like methods the amplitude of the full deviation is 
about 10~^, which coincides with the corresponding deviation of coy{x,Vx} shown in Fig. HI 
This indicates that, for the Gear-like methods, the errors in calculation come mainly from the 
evaluation of variances and covariances, that is, from the addition of the ERDs. On the other 
hand, for the Beeman-like method and particularly for the Euler-like method, the magnitudes of 
the full deviations are larger than those of the deviations of variances and covariances shown in 
Fig.m This implies that the errors introduced by addition of the ERDs might have been amplified 
during the calculation of the means. All in all, the dramatic differences in deviations in the x and 
z directions seen between different methods in Figs. [3l|4] and [5] have disappeared when ERDs are 
added. 
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FIG. 10: (Color online) Deviations of x and Vx (left), and z and Vz (right) from the corresponding analytical 



solutions i40h using different methods (from top to bottom: GL-5, GL-4, GL-3, BL and EL) in 200 time 



units, with 7 = 0.02, ujq = l/\/2, = 0.5 and = 0.02. Initial conditions are the same as in Fig. [T] 
Note the magnitude of the errors. 

From now on, we shall examine mostly the behavior of simulations by monitoring the quantity 

E{t)=v\t)+uy{t), (31) 



where v"^ = v1 + Vy + vl and 



+ y"^ + z^, which is defined to be proportional to the 



total energy of the Brownian oscillator at time t and, as such, it contains inaccuracies in both 
the position and velocity. Note that, because of coupling with the medium through damping and 
Brownian acceleration, this energy is no longer a conserved quantity. In order to examine the 
energy conservation performance of our numerical methods, we normalize E{t) by its "exact" 
counterpart Ee{t), which is simply obtained from Eq. (|3TI) by substituting the explicit analytical 
expressions for the position and velocity ll40ll . The resultant ratio E{t)/Ee{t) should be then a 
conserved quantity in a simulation with the expected value of unity. 



Figure \TT\ displays the relative deviation of energy from its exact counterpart i40|] . i.e., 
E{t)/Ee{t) — 1, calculated by using different methods for several damping rates, with At = 0.02 
and f2 = 0.5. To check the long time stability of these methods, a longer time scale of 1000 time 
units is adopted here. One sees that, for the Gear- like methods, the amplitude of deviation is quite 
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FIG. 11: Deviations of the total energy from its exact counterpart 14011 . i.e., E{t)/Ee{t) — 1, using different 
methods (from top to bottom: GL-5, GL-4, GL-3, BL and EL) for different damping rates 7 in 1000 time 
units, with = l/\/2, ^ = 0.5 and At = 0.02. Initial conditions are the same as in Fig. [T] Note the 
magnitude of the errors. 

similar to that of the position and velocity in Fig. \W\ Also, the amplitude does not change ap- 
preciably with the damping rate, although the frequency of the noise changes dramatically as the 
collision frequency, i.e., the damping rate, increases from 7 = 0.01 to 1.0. However, the situation 
for the Beeman-like and Euler-like methods is quite different. For the former, the amplitude of 
the deviation is approximately one order higher than that of the Gear-like methods at 7 = 0.01, 
but it decreases dramatically when 7 increases, and reaches almost the same level as that of the 
Gear-like methods at about 7 = 1.0. For the latter, the amplitude of the deviation is about 1 at 
7 = 0.01, and is therefore comparable to the value of energy itself, indicating that this method is 
not stable under these conditions. However, it also decreases with increasing damping rate. So, it 
is obvious that finite damping 7 actually stabilizes Beeman-like and Euler-like methods. This is 
not surprising at all, because the damping could also diminish errors inherited from a previous step 
during simulation. Indeed, previous studies 113 811 have demonstrated a possibility of stabilizing the 
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FIG. 12: (Color online) Largest deviation in total energy from its exact counterpart 114011 in 1000 time units 
versus At with cjq = 1/ a/2 and Q = 0.5 for different methods and different 7 values. Initial conditions are 
the same as those in Fig.[T] 



MD simulation by introducing a small, but finite damping. 

A more quantitative analysis of these results is shown in Fig. [121 in which the largest deviations 
in energy during 1000 time units are plotted versus the time step size At for different 7 values. 
One can see that the deviation of the Euler-like method is always the highest, but it decreases with 
increasing 7. Its magnitude suggests that this method should only be used for simulating systems 
with damping rate larger than 1.0 and with very small time steps. As for the Gear-like methods, 
judging from the magnitude of deviations and the apparent order of around 2.5, their accuracy has 
reached the limit determined by the accuracy of deviations in variances and covariances shown 
in Fig. [9l Therefore, their deviations in the energy remain very close to each other, and they are 
nearly independent of 7. The performance of the Beeman-like method is in between the Euler-like 
and Gear-like methods for very small 7, but it reaches the same level as the Gear-like methods for 
large damping rates, such as 7 = 1.0. 



VII. CONCLUSIONS 



We have presented several new algorithms for studying Brownian dynamics of charged particles 
in an external magnetic field. All these methods were tested by comparison with the available 
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analytical results for a three-dimensional, Brownian-harmonic-oscillator model in the presence of 



an external magnetic field 114011 . It was found that the Gear-like method generally has the best 
performance in terms of accuracy, long time stability, and energy drift in a wide range of damping 
rates, and especially in the low-damping limit. Therefore, the Gear-like method should be highly 
recommended when studying systems with very low-damping and/or when using the BD method 



as a thermostat 113711 in a MD simulation. The Beeman-like method can also cover a wide range 
of damping rates with reasonably good accuracy and with negligible energy drift. It should be 
recommended for simulating systems with intermediate damping rates. The Euler-like method, as 
can be expected, has the poorest performance and can be used with confidence only in simulating 
over-damped systems, such as colloidal suspensions and/or polymeric fluids. Further detailed tests 



based on applications to magnetized complex plasmas will be presented elsewhere ll39ll 



We note that, besides applications in plasma physics, our numerical method could be also of 



interest in numerical studies of some stochastic processes in statistical physics 



HQ 



4411 . since we have actually tested here the recently developed analytical model for a Brownian- 



harmonic-oscillator in the presence of a magnetic field. 
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